#ifndef CUDA_RAND_H
#define CUDA_RAND_H

struct RandGen{
	unsigned int z1, z2, z3, z4;  
	
	RandGen(){}

	RandGen(unsigned int _z1,
		unsigned int _z2,
		unsigned int _z3,
		unsigned int _z4){
	
		z1 = _z1; z2 = _z2; z3 = _z3; z4 = _z4;
	}

	__host__ __device__
	unsigned int TausStep(unsigned &z, int S1, int S2, int S3, unsigned M)  
	{  
  		unsigned b=(((z << S1) ^ z) >> S2);  
  		return z = (((z & M) << S3) ^ b);  
	}

	__host__ __device__
	unsigned int LCGStep(unsigned &z, unsigned A, unsigned C)  
	{  
  		return z=(A*z+C);  
	}    

	__host__ __device__ 
	unsigned int cudaRandInt(){
		return (TausStep(z1, 13, 19, 12, 4294967294UL) ^  // p1=2^31-1  
		TausStep(z2, 2, 25, 4, 4294967288UL) ^    // p2=2^30-1  
		TausStep(z3, 3, 11, 17, 4294967280UL) ^   // p3=2^28-1  
		LCGStep(z4, 1664525, 1013904223UL));        // p4=2^32  
	}

	__host__ __device__
	double cudaRandUni(){
		return 2.3283064365387e-10 * cudaRandInt();
	}
	
	__host__ __device__
	double cudaNormal(){
		double u0 = cudaRandUni();
		double u1 = cudaRandUni();
		double r = sqrt(log(u0)*(-2));
		double theta = 2*3.1415927*u1;
		return r*sin(theta);
	}

};

#endif
